home *** CD-ROM | disk | FTP | other *** search
/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / aminet / misc / sci / ephem_src_4_28.lha / cal_mjd.c < prev    next >
C/C++ Source or Header  |  1992-04-17  |  3KB  |  148 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "astro.h"
  4.  
  5. /* given a date in months, mn, days, dy, years, yr,
  6.  * return the modified Julian date (number of days elapsed since 1900 jan 0.5),
  7.  * *mjd.
  8.  */
  9. cal_mjd (mn, dy, yr, mjd)
  10. int mn, yr;
  11. double dy;
  12. double *mjd;
  13. {
  14.     int b, d, m, y;
  15.     long c;
  16.  
  17.     m = mn;
  18.     y = (yr < 0) ? yr + 1 : yr;
  19.     if (mn < 3) {
  20.         m += 12;
  21.         y -= 1;
  22.     }
  23.  
  24.     if (yr < 1582 || yr == 1582 && (mn < 10 || mn == 10 && dy < 15)) 
  25.         b = 0;
  26.     else {
  27.         int a;
  28.         a = y/100;
  29.         b = 2 - a + a/4;
  30.     }
  31.  
  32.     if (y < 0)
  33.         c = (long)((365.25*y) - 0.75) - 694025L;
  34.     else
  35.         c = (long)(365.25*y) - 694025L;
  36.  
  37.     d = 30.6001*(m+1);
  38.  
  39.     *mjd = b + c + d + dy - 0.5;
  40. }
  41.  
  42. /* given the modified Julian date (number of days elapsed since 1900 jan 0.5,),
  43.  * mjd, return the calendar date in months, *mn, days, *dy, and years, *yr.
  44.  */
  45. mjd_cal (mjd, mn, dy, yr)
  46. double mjd;
  47. int *mn, *yr;
  48. double *dy;
  49. {
  50.     double d, f;
  51.     double i, a, b, ce, g;
  52.  
  53.     d = mjd + 0.5;
  54.     i = floor(d);
  55.     f = d-i;
  56.     if (f == 1) {
  57.         f = 0;
  58.         i += 1;
  59.     }
  60.  
  61.     if (i > -115860.0) {
  62.         a = floor((i/36524.25)+.9983573)+14;
  63.         i += 1 + a - floor(a/4.0);
  64.     }
  65.  
  66.     b = floor((i/365.25)+.802601);
  67.     ce = i - floor((365.25*b)+.750001)+416;
  68.     g = floor(ce/30.6001);
  69.     *mn = g - 1;
  70.     *dy = ce - floor(30.6001*g)+f;
  71.     *yr = b + 1899;
  72.  
  73.     if (g > 13.5)
  74.         *mn = g - 13;
  75.     if (*mn < 2.5)
  76.         *yr = b + 1900;
  77.     if (*yr < 1)
  78.         *yr -= 1;
  79. }
  80.  
  81. /* given an mjd, set *dow to 0..6 according to which dayof the week it falls
  82.  * on (0=sunday) or set it to -1 if can't figure it out.
  83.  */
  84. mjd_dow (mjd, dow)
  85. double mjd;
  86. int *dow;
  87. {
  88.     /* cal_mjd() uses Gregorian dates on or after Oct 15, 1582.
  89.      * (Pope Gregory XIII dropped 10 days, Oct 5..14, and improved the leap-
  90.      * year algorithm). however, Great Britian and the colonies did not
  91.      * adopt it until Sept 14, 1752 (they dropped 11 days, Sept 3-13,
  92.      * due to additional accumulated error). leap years before 1752 thus
  93.      * can not easily be accounted for from the cal_mjd() number...
  94.      */
  95.     if (mjd < -53798.5) {
  96.         /* pre sept 14, 1752 too hard to correct */
  97.         *dow = -1;
  98.         return;
  99.     }
  100.     *dow = ((long)floor(mjd-.5) + 1) % 7;/* 1/1/1900 (mjd 0.5) is a Monday*/
  101.     if (*dow < 0)
  102.         *dow += 7;
  103. }
  104.  
  105. /* given a mjd, return the the number of days in the month.  */
  106. mjd_dpm (mjd, ndays)
  107. double mjd;
  108. int *ndays;
  109. {
  110.     static short dpm[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
  111.     int m, y;
  112.     double d;
  113.  
  114.     mjd_cal (mjd, &m, &d, &y);
  115.     *ndays = (m==2 && ((y%4==0 && y%100!=0)||y%400==0)) ? 29 : dpm[m-1];
  116. }
  117.  
  118.  
  119. /* given a mjd, return the year as a double. */
  120. mjd_year (mjd, yr)
  121. double mjd;
  122. double *yr;
  123. {
  124.     int m, y;
  125.     double d;
  126.     double e0, e1;    /* mjd of start of this year, start of next year */
  127.  
  128.     mjd_cal (mjd, &m, &d, &y);
  129.     if (y == -1) y = -2;
  130.     cal_mjd (1, 1.0, y, &e0);
  131.     cal_mjd (1, 1.0, y+1, &e1);
  132.     *yr = y + (mjd - e0)/(e1 - e0);
  133. }
  134.  
  135. /* given a decimal year, return mjd */
  136. year_mjd (y, mjd)
  137. double y;
  138. double *mjd;
  139. {
  140.     double e0, e1;    /* mjd of start of this year, start of next year */
  141.     int yf = floor (y);
  142.     if (yf == -1) yf = -2;
  143.  
  144.     cal_mjd (1, 1.0, yf, &e0);
  145.     cal_mjd (1, 1.0, yf+1, &e1);
  146.     *mjd = e0 + (y - yf)*(e1-e0);
  147. }
  148.